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We develop a mathematical model of the peripheral naive T cell population to study the 
change in human naiveT cell numbers from birth to adulthood, incorporating thymic output 
and the availability of interleukin-7 (IL-7). The model is formulated as three ordinary differ- 
ential equations: two describe T cell numbers, in a resting state and progressing through 
the cell cycle. The third is introduced to describe changes in IL-7 availability. Thymic output 
is a decreasing function of time, representative of the thymic atrophy observed in aging 
humans. Each T cell is assumed to possess two interleukin-7 receptor (IL-7R) signaling 
thresholds: a survival threshold and a second, higher, proliferation threshold. If the IL-7R 
signaling strength is below its survival threshold, a cell may undergo apoptosis. When the 
signaling strength is above the survival threshold, but below the proliferation threshold, 
the cell survives but does not divide. Signaling strength above the proliferation threshold 
enables entry into cell cycle. Assuming that individual cell thresholds are log-normally dis- 
tributed, we derive population-average rates for apoptosis and entry into cell cycle. We 
have analyzed the adiabatic change in homeostasis as thymic output decreases. With a 
parameter set representative of a healthy individual, the model predicts a unique equilib- 
rium number of T cells. In a parameter range representative of persistent viral or bacterial 
infection, where naiveT cell cycle progression is impaired, a decrease in thymic output 
may result in the collapse of the naiveT cell repertoire. 
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1. INTRODUCTION 

The number of naive T cells in the periphery is determined by a 
balance between cell loss (death or differentiation) and cell renewal 
due to cell division and thymic export (1,2). In humans, at least, 
the decline in thymic export occurs mainly in childhood, from 
about a year of age until 20 years of age, when the number of naive 
T cells is increasing (3). In adults, the decline in thymic export is 
much less pronounced but the number of naive T cells is, more 
or less, constant (4). Survival of the naive T cell population in the 
periphery depends on both common gamma chain cytokines and 
weak "tonic" signals induced by recognition of self-peptides by the 
T cell receptor (TCR) (5, 6). IL-7 is required for the homeostatic 
expansion of naive CD8 + and CD4 + T cells in lymphopenic hosts, 
while naive T cells disappear over a 1 -month period upon adoptive 
transfer into IL-7 deficient (IL-7 - ) hosts (7-9). 

Signals from recognition of self-peptides bound to major histo- 
compatibility complex (sp-MHC), and IL-7, promote cell survival. 
Naive T cell survival is impaired when removing access to one of 
these signals (10-14). Of interest are the mechanisms by which 
these signals are regulated, and that result in a stable number of 
naive T cells throughout the lifetimes of mice and humans. In this 
paper, we focus on IL-7 as a master regulator of naive T cell sur- 
vival (15). IL-7 is produced by stromal cells in tissues, including 
fibroblastic reticular cells, marginal reticular cells, and lymphatic 
endothelial cells (16). These cells produce very small amounts of 
IL-7 messenger RNA, consistent with IL-7 protein levels limiting 



T cell expansion. IL-7 is a heparin-sulfate binding protein, and as 
such, it will bind extra-cellular matrix surrounding stromal cells. 
Thus, the interaction between naive T cells and stroma controls 
their homeostasis (17). Recognition of higher affinity, non-self- 
peptides by the T cell receptor induces naive T cells to undergo 
an alternative, IL-7 independent, survival program dependent on 
IL-2 (18). 

Naive CD8 + T cell responses depend on the amount of IL-7 cells 
are exposed to (19). At low IL-7 concentrations (<10 -2 ngml -1 ), 
cell viability was impaired; at higher concentrations ( > 1 ng ml - 1 ) 
cells were observed to proliferate in response to IL-7. This dif- 
ference might arise from changes in the strength of the IL-7R 
induced signal the cell receives. For an individual cell, IL-7R 
induced signaling must be greater than some threshold to prevent 
the accumulation of pro-apoptotic proteins. Similarly, IL-7R sig- 
naling must be greater than a second, higher, threshold to induce 
cell division. Heterogeneity at the single cell level in IL-7 signaling 
thresholds (a property reported to depend on expression of IL- 
7R), resulted in differential survival and division (19). Although 
these observations are based on two different CD8 + T cell receptor 
transgenic mice, it is assumed that the key principles regarding T 
cell survival will be found in the repertoire of naive CD4 + and 
CD8+ T cells. 

We introduce a deterministic mathematical model of the naive 
T cell population to study the change in human naive T cell 
numbers from birth to adulthood. We will assume cell survival 
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depends on the availability of IL-7. We do not include availability 
of sp-MHC as a variable within the model, but assume sp-MHC 
availability is sufficient to allow cell survival and proliferation, 
in conjunction with sufficient IL-7 stimulus. We also make the 
approximation that heterogeneity is constant with changes in age. 
For a mathematical study of the impact sp-MHC availability has 
on clonal diversity, the reader is referred to Stirk et al. (20, 21). 
Our model is a mathematical description of the homeostasis of 
the naive T cell repertoire, but does not consider stimulation by 
foreign antigens. 

2. MATERIALS AND METHODS 

2.1. A MATHEMATICAL DESCRIPTION OF THE SIZE OF THE 
PERIPHERAL NAIVE T CELL POPULATION 

Stochastic processes provide a method of treating each cell as a dis- 
tinct, countable object, and permit a more realistic model than a 
deterministic characterization. Fluctuations in the number of cells 
can be considered but, in a non-linear stochastic model, approxi- 
mations are often made to facilitate the analysis. In the linear noise 

approximation (22), for example, fluctuations are assumed to be 
l 

of order £22 for a system of size Q. The human peripheral T cell 
compartment is estimated to contain of the order of 10 11 T cells 
(3). Letting the system size be the average number of naive T cells 
in humans, we find 0(Q) = 10 11 cells, and correspondingly, fluc- 
tuations are expected to be typically 10 5 — 10 6 cells in magnitude. 
That is, we expect fluctuations of approximately 0.001% in the size 
of the human naive T cell pool due to stochasticity in the per cell 
division and death rates. Based on these considerations, adopting 
a deterministic approach to describe the total human peripheral 
naive T cell population is reasonable. 

We assume peripheral naive T cells are either in a resting state, 
or proceeding through the cell cycle. The deterministic variables 
R(t) and C(t) are introduced to model the total number of T cells 
in the resting and cycling states, respectively. The variable I(t) is 
introduced to model the concentration of IL-7. The deterministic 
approach we take does not consider any notion of space. Indeed, 
this approach is tantamount to assuming the resource, IL-7, is 
shared equally amongst all cells. Competition for the resource is 
introduced only so far as each cell acts to reduce the global con- 
centration of the resource. Resting cells may receive a signal which 
induces them to proceed through one round of division. Upon 
completion of the cell cycle, a cycling cell produces two daughter 
cells in the resting compartment. Resting cells are assumed to die 
if the IL-7 induced survival signal is insufficient; cells may also 
die during cell cycle. The input of cells from the thymus into the 
resting compartment, in keeping with observations in humans, is a 
decreasing function of time (23, 24). Production of IL-7 is related 
to the size of the lymphatic system architecture, which we estimate 
from the body mass of an individual. In the absence of T cells, IL-7 
is assumed to be degraded and/or consumed by other cell types at 
a constant rate. Upon signal induction through the IL-7 receptor, 
IL-7 is assumed to be consumed by the T cell. A diagrammatic 
representation of the model is given in Figure 1. 

2.2. IL-7 SIGNALING AND HETEROGENEITY IN IL-7 THRESHOLDS 

In the model, IL-7 signaling is assumed to be uniform across 
the population. Yet, we introduce heterogeneity in the signaling 
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FIGURE 1 | Diagrammatic illustration of the deterministic model. T cells 
leaving the thymus enter the resting naive peripheral pool. Cells in either a 
resting or cycling state may die. The rate of death from the resting state 
depends on the availability of the resource (IL-7), whereas the death rate for 
cycling cells is constant. Resting cells enter the cell cycle at a rate that 
depends on the availability of the resource (IL-7). Cycling cells produce two 
daughter cells in the resting state upon completion of the cell cycle. 



thresholds for survival and proliferation. Let S(t) be the aver- 
age signaling strength across the naive T cell population. Each 
cell experiences the same strength of signaling for a given con- 
centration of IL-7, I(t). We relate the internal signaling to the 
concentration of IL-7 by the equation 



S(t) 



600/ ( t) 
0.025 + 1 (t) 



(1) 



The functional form and constants of this relationship are 
derived from the study of IL-7 receptor dynamics summarized in 
the Appendix. We assume each individual cell possesses a unique 
pair of thresholds for survival and proliferation. Furthermore, we 
assume, in the continuous limit, that these thresholds are distrib- 
uted log-normally across the entire population of cells. The use of 
log-normal distributions guarantees, first of all, that all signaling 
thresholds are positive real numbers. Secondly, the log-normal dis- 
tribution ensures that no cell can survive or divide independently 
of IL-7. 

Let the random variable @ 5 represent the survival threshold, 
and let ©p represent the proliferation threshold. We write 

0 S ~ logA/" (k)g£ s , , O p ~ logA/" (k)g£ p , , 



a e . 



(2) 



The respective probability density functions are 
p@ x (0) = -?L= exp \-(a (\og0 - log^)) 2 ] , x = s, p. (3) 

Estimates for 6 S and 0 p are obtained from the model summa- 
rized in the Appendix. Lauffenburger et al. found a significant 
change in cell viability for both OT-1 and F-5 T cells at around 
10 _25 ngml _1 IL-7 (19). Proliferation of OT-1 cells occurred 
above lngml -1 IL-7. We estimate the equilibrium signaling at 
these concentrations as 60 and 600 units, respectively, setting 
0 S = 60 and Op = 600. Modeling heterogeneity in IL-7 responses 
by assuming heterogeneity in the IL-7 signaling thresholds, allows 
us to avoid modeling the naive population using either: (i) a PDE 
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approach, where heterogeneity is continuous across the popula- 
tion of cells, or (ii) describing each subset of cells sharing common 
thresholds with its unique governing set of ODEs. However, these 
approaches present an obvious avenue for further research beyond 
the scope of this paper. 

2.2. 1. Death rate of resting cells 

Suppose each T cell in the naive population is a distinct mem- 
ber possessing a unique signaling threshold for survival. As in the 
previous section, we assume these signaling thresholds are dis- 
tributed log-normally (in the continuous limit). We suppose the 
death rate of an individual cell is Boolean, in the sense that if 
the global signaling strength is greater than the cell's individual 
survival threshold, then the cell can survive indefinitely. Similarly, 
if the signaling strength is below the survival threshold, the cell 
will undergo apoptosis at a rate /jLr. The death rate for cell i with 
survival threshold 0^ , is given by 



fs(s(t),0®) 



PR 
0 



if S (t) < 6> s (i) 
if S (t) > 0^ , 



(4) 



In the continuous limit (assuming signaling thresholds are dis- 
tributed log-normally), the average death rate for the population 
of naive T cells is given by 

poo poo 

Afi (S (t)) = / fs(S (t) , 6) p @s (0) d9 = / n R p @s (9) 69, 
Jo JSU) 



-Wi[l-erf(a(logS(f)-log0,))] , 



(5) 



where p@ s (6) is the probability density function of the random 
variable 0 S , defined by equation (3) with x = s. 

2.2.2. Rate of entry into cell cycle 

Analogous to Section 2.2.1, we assume each T cell in the naive 
population is a distinct member possessing a unique signaling 
threshold for proliferation. We let the individual rate of entry into 
cell cycle be given by 



f p (s(t), 0®) 



0 ifS(0< 
P if S(t) > 0, 



a(0 
(i) 



(6) 



Assume, in the continuous limit, the signaling threshold for 
entry into the cell cycle is represented by the random variable 0p, 
defined in equation (2), with probability density function p@ p (0) 
(equation (3), x = p). The average rate of entry into cell cycle is 
given by 

poo pS{t) 

P (S (0) = / f P (S (0 , 0) p @p (0) dO = / pp @p (0) dO , 
Jo Jo 



-p[l + erf (a(logS(t)-log^))] . 



(7) 



2.2.3. Cell cycle progression 

Cycling cells take on average A -1 days to complete the cell cycle. 
After a cell divides, both daughter cells are produced in the resting 
state and require a second signal before they can progress through 
another round of cell division. Cell cycle may be interrupted 
resulting in the death of the cell. Such death events occur at a 
rate fie- 

2.2.4. Thymic export 

We assume thymic output to be a decreasing function of time. In 
particular, we use the functional form given by Bains et al. (3). Let 
us introduce the thymic output function, v(t), as follows 



v (f) = 2.32 x 10 8 exp (-1.1 x 10 4 f) 

+ 1.15 x 10 8 exp (-1.6 x 10 7 f 2 ) 



(8) 



where t corresponds to the age of the individual, measured in 
days. A plot of this function is shown in the left panel of Figure 2. 
The function was chosen by Bains et al. to describe the rate of 
thymic export of CD4 + T cells. We use the same function to 
describe the export rate of all naive T cells (CD4 + or CD8 + T cells) . 
This approximation is justified since we require the absolute cell 
count to roughly approximate the cell count observed in humans 
(indeed, such an observation is likely subject to large differences). 
Of interest later in the paper is the relative variation of cell num- 
bers with different choices of parameter values. For our purposes, 
the important feature of v(t) is that it is a decreasing function of 
time. 

2.2.5. Internalization of 11-7 

We use the model summarized in the Appendix to estimate the rate 
of IL-7 internalization. The total number of IL-7 molecules inter- 
nalized by a single T cell in 1 day, exposed to IL-7 at concentration 
I ngml -1 , is described by the function 



3 + 



x 10 3 cell -1 day -1 . (9) 



6.71 (f) 

2 + 1(0 1(0 + 3 x 10 

It is reported IL-7 has a molecular mass of 17kDa 
(^2.8 x 10 -11 ng) (25). Based on this, we define the per cell rate 
of IL-7 internalization to be 



2/(0 
2 + 1(0 



3 + 



1(0 + 3 x 10" 



x 10 7 ng day cell 



In order to convert the rate of change of mass to the rate of 
change of concentration, we must choose a volume for the system. 
Naive T cells are typically found in the lymph nodes, spleen, and 
gut of the human body. We shall make the rough estimation that 
the total volume is of the order of 1 1. This implies the rate of IL-7 
internalization by all naive T cells in the population is given by 

21(0 / 5 \ 

y (I (0) R (0 = — 3 + 

Y 2 + 7(0 V 1(f) + 3xl0- 2 / 

x 10 -10 £ (0 ng ml" 1 day -1 , (10) 

where R(t) is the number of resting naive T cells. The IL-7 
internalization rate, y(I(t)), is shown in the middle panel of 
Figure 2. 
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FIGURE 2 | Left panel: rate of thymic export. Plot of equation (8). Middle panel: IL-7 internalization rate. Plot of the function y{l{t)). Right panel: plot of the 
function m{s) (equation (11)), a model for the average body mass of male individuals. 



2.2.6. Production of IL-7 

We assume the rate of IL-7 production is proportional (with pro- 
portionality constant ft) to the body mass of an individual. To 
estimate average body mass we use the model given by Burmaster 
and Crouch (26). The explicit relationship between mass and age 
is given by the function m(s) as follows 

m(s) = exp [4.1 + 1.4 x 10" 2 s - 1.5 x 10~ 4 s 2 

-2.0exp[-s(0.15 - 1.4 x 10" 2 s + 9.8 x 10" 4 s 2 )]] 



(11) 



where s is age measured in years. A plot of this function is given 
in the right panel of Figure 2. The rate of IL-7 production is given 
by the function 



P(t) = 0m (365f), 
where t denotes age as measured in days. 



2.2. 7. Intra- cellular degradation of IL-7 

We assume IL-7 is degraded and internalized by other cell types at 
a constant rate. We let the degradation rate of IL-7 be denoted by 
the parameter 8. 

2.2.8. Deterministic mathematical model 

From the above assumptions, the system of differential equations 
governing the behavior of the naive T cells (resting and cycling) 
and the concentration of IL-7 is given by 



dl(t) 
dt 

dR(t) 
dt 

dC(t) 
dt 



= P(t)-y(I(t))R(t)-8I(t) , (13) 
= v (t) - [p (S (0) + [Ir (s (0)] R (t) + 2XC (t) , 



= p(S(t))R(t)-Qi c + k)C(t) 



(14) 
(15) 



3. RESULTS 

3.1. PARAMETER ESTIMATES 

The function describing the rate of internalization of IL-7, equa- 
tion (10), and the signaling relation, equation (1), were both esti- 
mated from our studies of IL-7 receptor dynamics in naive T cells 
(summarized in the Appendix). In the same study we have esti- 
mated the average signaling thresholds to be, respectively, 0 S = 60 



and^ 



: 600. Our estimate for 0 P was found to be close (6 P ~ 585, 



This system is subject to the initial conditions Iq, #o> and Cq. 



which we have rounded to 600) to the limit of the signaling func- 
tion equation (l),asI(f)^+oo. These estimates imply that, on 
average, approximately half of the naive T cell population does 
not proliferate in response to excess IL-7 when receptor dynamics 
is in equilibrium. Such a finding is consistent with the observa- 
tions made by Lauffenburger et al. in the experiments described 
in Ref. (19). In these experiments, F-5 T cells did not proliferate 
(whereas OT-1 cells did proliferate), even in excess amounts of IL- 
7. Estimates for the IL-7 average signaling thresholds, 0 S and Op, 
have been based on OT- 1 cells. However, we note that changes to 
these thresholds only result in quantitative differences, provided 
#s < @p- 

We assume cycling naive T cells take 12 h to complete the 
cell cycle and produce two daughter cells (A -1 =0.5 day). Rest- 
ing naive T cells are assumed to die after 2 days of IL-7 starva- 
tion (/x^ 1 = 2 day). We, furthermore, set the rate of entry into 
cell cycle to be p = 5day -1 . Cycling cells are assumed to die 
at rate [ic— 1 day -1 in healthy individuals, whereas we choose 
Ijl c = 3 day -1 (>A) to be representative of cell cycle impairment. 
These parameters have been estimated from the literature. We note 
that the model behavior we discuss in the following sections was 
found to be robust to changes in these parameters. Robustness 
was concluded since a 10-fold change in each or any combination 
thereof, of these four parameters did not change the qualitative 
behavior of the model, provided the relation A>/xc(orA.</xc) 
was maintained. We examine the changes in model behavior aris- 
ing from altering the relative values of A and /xc in the following 
sections. 

We choose 8 such that 8I(t) is similar in magnitude to 
y(I(t))R(t), when R(t)^ I0 l \ and I(t) « 10" 2 ngml" 1 . The 
proportionality constant /3 is chosen such that we observe 
O(l0 n ) naive T cells in equilibrium at 20 years of age (for 
I(t) ~ 10 -2 ng ml -1 ). We had the least inclination when choosing 
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Table 1 | Parameter choices for the mathematical model. 



Parameter 


Value 


Units 


0 S 


60 


Signaling units 




600 


Signaling units 


a 


2 


(Log signaling units) -1 




0.2 


10" 12 m|- 1 day" 1 


8 


500 


day -1 


X 


2 


day -1 


P 


5 


day -1 




0.5 


day -1 




1 


day -1 



The problem of finding adiabatic solutions can then be reduced 
to finding solutions, I(t), to the one-dimensional equation 



the parameter a. This parameter describes the spread in the indi- 
vidual signaling thresholds across the naive T cell population. We 
choose a = 2, however this choice has no justification from the 
literature. The parameter set is summarized in Table 1. 

3.2. THERE EXISTS A UNIQUE AND STEADY STATE WHEN X > ti c 

Let us suppose changes in thymic output and IL-7 production 
occur in time scales slower than those of the changes in the num- 
ber of naive T cells. Under this assumption, we look for adiabatic 
solutions of the system as follows: 

o = p(t)-y(i(t))R(t)-8i(t) , (16) 

0 = v (f) - (p (S (0) + Plr (S (t))) R (t) + 2XC (f) , (17) 
0 = p (s (0) R (t) - (ji c + A) C (f) . (18) 

For the parameter set studied, the relative error between this 
solution and the exact solution is within 3% for the resting naive 
T cell population and the concentration of IL-7, and within 14% 
for the cycling T cell population (see Figure 3). 

All numerical results presented in the paper have been obtained 
with a Python code 1 : differential equations (13), (14), and (15) 
have been solved using a fourth-order Runge-Kutta scheme. 
Quasi- stationary solutions were found using the scipy. optimize 
package. Bifurcation plots were computed using a bisection 
scheme to search for multiple solutions of equation (21) in the 
interval [0, 1 ] . Corresponding T cell numbers were calculated using 
equations (19) and (20). 

Notice that the adiabatic solution for both cell types is uniquely 
defined for a given value of cytokine concentration, I(f)> namely 



R(t) = 



C(t) = 



P(t)-8I(t) 

y(ht)) 

p (s (*)) p ( t ) _ si (t) 
^c + x y(l (o) 



(19) 



(20) 



v (t)y (1(0) 
P(t)-8I(t) 



= Ar (§(*)) 



+ 1 



2k 



A. + Mc 



p(S(0) . (21) 



1 Python code available upon request. 



By construction, y (1(f)) is a monotonically increasing function 
of 1(f). Furthermore, the existence of positive adiabatic solutions 
requires /3(t) > 8I(t). Therefore, for a fixed time f, the left-hand 
side of equation (21) is an increasing function of I(t). For X > fie, 
the right-hand side of equation (21) is a monotonically decreasing 
function of 1(f). Lastly, y (1(f)) = 0 for 1(f) = 0, and the limit as 
1(f) —> 0 of the right-hand side is equal to fiR. It follows that the 
intersection of the left and right sides must be unique and positive. 
We deduce that for X > there exists a unique adiabatic solution 
to the system (see left plot of Figure 4). This solution is stable 
for the parameter set given in Table 1. We have also numerically 
explored parameter space, but have not found a parameter set for 
which this solution is unstable. In Figure 5 we present numeri- 
cal solutions to equations (13-15), computed with a fourth-order 
Runge-Kutta method implemented in Python. Initial conditions 
were chosen to be the adiabatic solutions of equations (16-18) 
at f = 0. 

3.3. THERE EXIST TWO STEADY STATES WHEN X < [i c 

Let us now consider the case X < fie- The right-hand side of equa- 
tion (21) is no longer a decreasing function of 1(f). We find 
there may exist up to three solutions 1(f), two of which may 
be stable simultaneously, whilst the third is unstable (see right 
plot of Figure 4) . Further examination (by numerically finding 
all solutions using the bisection method) of the model reveals a 
saddle-node bifurcation as thymic output changes with age (see 
left panel, Figure 6). We assume that individuals with (around) 
10 11 naive T cells are healthy, in the sense that they have a suf- 
ficient number of T cells to provide protection against immune 
challenges. When a cell in cycle is more likely to die than to 
produce two daughter cells, we define cell cycle progression to 
be impaired. Mathematically, this corresponds to X < fie- When 
cell cycle progression is impaired, the model predicts an indi- 
vidual may possess a healthy number of T cells, thereby being 
immuno- competent, up until the age at which the model bifur- 
cates. For the parameter set we have investigated, this bifurca- 
tion is inevitable given the estimated decline in thymic output 
established by Bains et al. (3). Indeed, for a given parameter 
set, from the known rate at which thymic output declines, one 
can estimate the age at which the model bifurcates. The bifur- 
cation results in a decrease in the naive T cell population size 
of approximately two orders of magnitude, that is, following the 
bifurcation we expect roughly 99 out of every 100 naive T cells to 
be lost. 

4. DISCUSSION 

In a healthy individual, it is reasonable to expect that naive T cells 
entering the cell cycle are more likely to complete division and pro- 
duce two daughter cells, than to die during the division process. 
Therefore, we suppose the parameter relation X < \xq represents a 
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Time (years) Time (years) 



FIGURE 3 | Relative error between the adiabatic solution and the exact solution, where initial conditions were chosen to be equal to the adiabatic 
solution at time t = 0. 
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FIGURE 4 | Left panel: for k> fi c (t = 25 years), the right-hand side 
of equation (21) is a decreasing function of /(f) (red line). There 
exists a unique, asymptotically stable solution /(f) found at the 
intersection of the solid black and red lines. Right panel: there exist 



three intersections between the red and solid black lines when A < /x c 
for t = 25 years. We require f}(t) > y(l(t)) for existence of stable 
solutions, therefore we neglect all intersections with the dashed 
black line. 
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FIGURE 5 | Numerical solutions to equations (13-15) for X > fi c . Initial 
conditions are chosen to be the adiabatic solutions of equations (16-18) at 
f = 0. We observe a large increase in the naive T cell population from birth to 



adulthood. During this period we see a 10 7 -fold increase in the number of 
cycling cells, resulting from increased IL-7 availability and reduced thymic 
output. IL-7 availability increases continually as the individual ages. 



healthy individual. We have shown, under this hypothesis, that the 
model allows a single asymptotically stable adiabatic solution. The 
total number of naive T cells per kilogram of body mass was found 



to increase over the first 18 years of life and decrease thereafter. The 
decline in this ratio for adults is seemingly a consequence of the 
reduction in thymic output. 
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FIGURE 6 | Left panel: bifurcation diagram for varying the parameter v(t) for X < fi c . Right panel: bifurcation plot as a function of varying the rate of IL-7 
production. Thymic export is fixed at v(60 years). 



The number of cycling naive T cells in adiabatic conditions was 
found to increase from practically 0 (approximately 0.75 cells) 
cells at birth to 10 7 cells by adulthood. Whilst the number of 
cycling cells increased thereafter, the increase was at a markedly 
slower rate. The increase in the number of cycling cells is probably 
due to the decline in thymic output, wherein competition for IL-7 
decreases. Despite the increase in the cycling population, the num- 
ber of cells in cycle is several orders of magnitude smaller than the 
resting population. For adults, this proportion is approximately 
[0.01, 0.04]% of the total naive population. 

When cell cycle is impaired (X < /xc), the model exhibits a 
saddle-node bifurcation as thymic output declines. For the para- 
meter set investigated, the two adiabatic solutions are generally 
separated by two orders of magnitude. This feature motivates the 
following theoretical scenario: suppose a healthy individual expe- 
riences some event which causes naive cell cycle progression to 
become impaired at 18 years of age. The model predicts the naive T 
cell repertoire will experience a dramatic reduction in cell numbers 
at roughly 33 years of age. See Figure 7 for the full solution of the 
model illustrating this scenario. From an immunological perspec- 
tive, this scenario is interesting because the noticeable effect of the 
event (the dramatic loss of naive T cells) occurs roughly 15 years 
later than the event itself (naive cell cycle progression to become 
impaired at 18 years of age). Indeed, it is the reduction in thymic 
output that triggers the loss rather than the event itself. If the 
thymus was not atrophic, the loss in naive T cells would not occur. 

Suppose now we fix thymic output at a constant rate and let 
the cell cycle be impaired. The model undergoes a saddle-node 
bifurcation as the rate of IL-7 production changes. More specif- 
ically, the model bifurcates, resulting in a decrease in the total 
number of T cells, as a consequence of increasing IL-7 produc- 
tion. Intuitively, this can be understood as follows: for increased 
IL-7 production, the rate of entry into cell cycle is enhanced, 
however, since cells are more likely to die in cell cycle than to 
produce two daughter cells, the enhanced rate of entry into cell 
cycle actually serves to decrease the total amount of naive T cells. 
The bifurcation diagram for this scenario is shown in the right 
panel of Figure 6. We found the critical value of this bifurcation 
decreases with thymic output. In the limiting case, corresponding 
to thymic output at 60 years of age, we found this critical value 
to be approximately 9.1ngml -1 . Consider again the theoretical 



scenario in which a healthy individual undergoes an event result- 
ing in cell cycle impairment at 18 years of age. Suppose now at 
age 25 the same individual undergoes some treatment to limit IL- 
7 production to 8ngml _1 , corresponding to thymic production 
at approximately 1 1 years of age. The T cell count is reduced by 
approximately 60%. This reduction is a significant improvement 
on the 99% T cell loss in the untreated individual at 33 years of 
age. For this theoretical scenario, the model predicts limiting IL- 
7 availability will partially avoid the dramatic T cell loss arising 
from reducing thymic export when the cell cycle is impaired. See 
Figure 8 for the full model solution in the treated individual. 

In the model presented here we have neglected the fact naive 
T cells become activated in response to recognition of ligand spe- 
cific to their unique TCR. Research by Koenen et al. has shown T 
cell survival is IL-7 independent following T cell activation (18). 
Suppose now we include a term in the governing ODEs to rep- 
resent differentiation of naive T cells into cells with a different 
phenotype, such as activated T cells. Assuming no reversion back 
to the naive phenotype, such a term would appear in the model 
as a loss term equivalent to the death rate /x#(S(f))- The simplest 
approach to including differentiation (due to activation) would 
be to assume naive T cells differentiate into activated T cells at a 
constant rate proportional to the rate of antigenic challenge. In 
this case, we would replace the term /x#(S(f)) by /x#(S(f)) + /X£>, 
where /X£> is constant. Consider again the red curves in Figure 4. 
The differentiation term will cause a translation in the red curve 
of length up the vertical axis. For A > /xc> there still exists a 
unique solution, however there will be a quantitative change in 
comparison to the case when /X£> = 0. When X < \ic, there still 
exists the possibility of a saddle-node bifurcation for a general 
parameter set. For the parameter set we have investigated, there 
exists a maximum value /x^, for a given time t, such that we can 
find more than one solution. For differentiation rates fio > /x^, 
we can only find a unique solution corresponding to the stable 
adiabatic one, in which we have reduced T cell numbers. 

In this paper we have developed and analyzed a determinis- 
tic mathematical model of a population of naive T cells, whose 
survival depends on the availability of the cytokine IL-7. We have 
shown this model predicts a stable population of cells when cell 
cycle progression is healthy. More interestingly, when cell cycle 
progression is impaired, our results indicate declining thymic 
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FIGURE 7 | Theoretical scenario in which we set fi c = 3 at 18 years of age. At roughly 33 years of age the model bifurcates resulting in a dramatic loss of 
naive T cells. 
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FIGURE 8 | At age 18 the individual experiences some event resulting in 
impaired cell cycle progression. At age 25 the individual undergoes 
treatment to limit IL-7 production to levels comparable to those at 11 years of 



age. Whilst theT cell count is reduced, this attrition is a significant 
improvement on the reduction observed when the individual has not received 
treatment (see Figure 7). 



output may result in a dramatic loss in the number of naive T 
cells. Furthermore, we have been able to establish that limiting IL- 
7 production partially rescues this decline. However, our study has 
been restricted to naive T cells and we have not taken into account 
the accumulation of memory T cells as an individual ages. Mem- 
ory T cells are generated in response to antigen or homeostatic 
cytokines (27), and during repeated homeostasis-driven divisions 
of naive T cells (28). Previous studies have shown the percent- 
age of memory T cells increases with age, yet the percentage of 
naive T cells decreases with age (29). CD4 + memory T cells also 
require IL-7 for survival and proliferation in the periphery (30). In 
the case of CD8 + memory T cells, IL-15 is largely responsible for 
their peripheral survival, but IL-7 may also be required. It is, then, 
reasonable to assume that memory T cells compete for the IL-7 
required for naive T cell survival. At least, that is, for those memory 
T cells which access IL-7 in the same tissues as naive T cells, such 
as the lymph nodes. As mentioned before, memory T cells can be 
derived from repeated divisions of naive T cells (3 1 ). In our model, 
the naive T cells that will acquire a memory-like phenotype first 
are those with low IL-7 survival and proliferation thresholds. We 
would expect, then, that over time the ratio of naive to memory 



T cells would decrease. Furthermore, the naive population would 
lose those cells with lowest survival and division thresholds. Given 
the distribution of signaling thresholds in our model, we would 
expect to see a shift to the right for the distribution of both survival 
and signaling thresholds. This shift, which implies that on average 
naive T cells require a higher concentration of IL-7 to enter cell 
cycle with age, together with the additional competition from an 
increasing memory population, is expected to cause a decrease in 
the total number of naive T cells and the percentage of naive T 
cells in cell cycle. These considerations might help explain the dis- 
crepancy between our model and data showing a reduction in the 
percentage of naive T cells expressing Ki67 + during childhood 2 . 
There is support for the fact that childhood might be the phase in 
which most memory T cells are acquired (3). Indeed, it might be 
that the model presented here, better describes the total number 
of T cells that require IL-7, naive and memory, rather than just 
the naive T cell pool alone. This, however, would require a more 
detailed analysis, out of the scope of this paper. 



Our current model suggests the percentage of naive T cells in cell cycle is increasing. 
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Many different receptor-mediated signals are integrated by T 
cells in their micro -environment, such as cytokines, adhesion 
molecules, T cell receptors, and co-receptors (CD28, CTLA-4). 
However, the survival of naive T cells, the focus of this paper, has 
been shown to depend on IL-7 and low level TCR stimulation 
(8, 9, 15, 32). Other cytokines, such as IL-15, play a significant 
role in the homeostasis of memory CD8 + T cells (33). Given the 
experimental support for the hypothesis that IL-7 is critical for the 
homeostatic proliferation and the survival of naive T cells, in this 
first model we have neglected other signals (7). 

Previous mathematical models of naive T cell homeostasis have 
focused on the relative contribution of thymic export and cell divi- 
sion in the periphery (3, 34, 35). These models conclude that in 
humans, thymic export makes an important contribution to the 
size of the naive T cell population in early life (< 20 years of age), 
whereas later in life the number of naive T cells is maintained 
by homeostatic proliferation in the periphery. Some recent stud- 
ies have measured the relative contribution of thymic export by 
examining the average number of TRECs in naive T cells (3, 35, 
36). The use of mathematical models has allowed these groups to 
infer the relative (young versus old) kinetics for naive T cells, based 
on experimental estimates of the total number of naive T cells and 
recent thymic emigrants (3, 34, 35). In the model introduced here, 
we have aimed to provide a mechanistic perspective by investi- 
gating at the molecular and cellular levels, the role IL-7 plays in 
regulating the homeostasis of naive T cells (37). The increase in the 
proportion of cycling cells in our model is in agreement with pre- 
vious experimental studies (38). Our study suggests the increase 
in peripheral division rates can be explained by the availability of 
IL-7, which is a consequence of a combined effect: (i) an increased 
net IL-7 production as an individual ages, and (ii) a reduction of 
recently exported thymocytes competing for this resource. 
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APPENDIX 

A.1 . A MODEL OF IL-7R DYNAMICS 

In this section we present a summary of a stochastic model in 
which we consider the number of IL-7 receptors on the surface of 
naive T cells. This study is used to derive equations (1) and (9) in 
the main text. The model has been formulated as a continuous time 
Markov process. We consider the number of receptors on a single 
naive T cell. For our purposes, we present the mean field approx- 
imation to the model with which we estimate the parameters. We 
introduce the following four variables: 

• m\(t) - the total number of IL-7 receptors on the surface of a 
naive T cell, 

• mi{t) - the number of unbound internalized receptors, 

• m$(t)- the number of bound internalized receptors, 

• m±{t)- the quantity of IL-7 induced signal. 

We assume, in equilibrium, the fraction of surface receptors 
bound to IL-7 is given by ft. The system of ODEs denoting the 
mean field approximation to the stochastic model is given by: 



drai (t) 
dt 



+ £t/ m 2 (t) + % B m 3 (t) 



dm 2 (0 
dt 

d/773 (t) 
dt 

dm4 (t) 
dt 



- [<xu(l -fi)+o B fi\ mi (t) , 
: cru (l - fi) mi it) - + Su) m 2 (t) 

■ CTBfl m (t) - (t; B + f B ) m 3 (f) , 

: (p m (t) - x m (t) . 



(Al) 
(A2) 

(A3) 

(A4) 



A.2. REDUCED MODEL IN THE CASE WHEN / = 0 

Consider a T cell in an IL-7 free medium with initial conditions 
such that the IL-7 induced signaling vanishes and the number of 
IL-7:IL-7R internal complexes is zero. Then the mean field model 
can be reduced to the following set of ODEs: 



drai (t) 

dt 
dm 2 (t) 

dt 



= 0 + §t/ m 2 (t) - o\j mi (t) , 
= cru m (t) - (£[/ + m 2 (t) . 



This reduced system is governed by four parameters 0, § u> & u> 
and £ u, and possesses the following stable steady state: 

* 0 (gu + ?u) 



?L7 



We assume in the reduced model 10% of the total number 
of receptors are internalized in equilibrium. We set m* = 9m\. 
Based on the measurements of Singer et al. (39), we shall assume 



4 x 10 4 receptors in total when the reduced model is in steady 
state. Therefore, we let 



»(§[/ + £[/) a . in 4 

= 3.6 x 10 , 



= 4x10'. 



(A5) 
(A6) 



In Ref. (40), cells were cultured with the translation inhibitor 
cycloheximide (CHX) to prevent transcription of the IL-7 recep- 
tor. Total expression of the IL-7 receptor was measured over several 
time points, from which the authors estimate the half-life of the 
receptor in an unstimulated cell to be approximately 24 h. 

In the reduced model, all receptors are guaranteed to be 
degraded in a finite amount of time. The expected time for a 
receptor, that is initially on the cell surface, to be degraded in the 
lysosome is given by 



?i = 



Assuming exponential decay, the half-life for a receptor to 
undergo lysosomal degradation, starting on the cell surface, is then 
given by 

ti = log 2. 

2 GUSU 



Thus, we can write 



W2 = 24 h. 



(A7) 



Combining equations (A5), (A6), and (A7) we find 

ft/ « 0.29 h" 1 , 
0 ~ 1.2 x 10 3 receptors h _1 , 
$u + 0.29 h" 1 ^9a v . 

The value of § u relative to f u effectively dictates the ratio of 
receptors which are degraded to those recycled back to the cell 
surface. We assume the system has evolved to minimize waste of 
functional proteins and tentatively let § u > f tj- That is, we assume 
that a greater fraction of receptors are recycled back to the surface 
of the cell. We somewhat arbitrarily set 



= 1 h" 



ojj « 0.14 h" 



A3. RECEPTOR-LIGAND KINETICS 

Suppose the number of surface receptors is constant and denoted 
by Rr. Let us also assume the extra- cellular concentration of IL-7 is 
constant and denoted by I. Define Rs(t) to be the number of IL-7 
receptors bound to IL-7. Note that we assume the time to recruit 
the common gamma chain, y c , is negligible. Then, we can describe 
changes in the number of bound complexes by the following ODE 



dR B (t) 
dt 



k+[R T -R B (t)]I-k-R B (t) 
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FIGURE A1 | Diagrammatic representation of the transition probabilities of the stochastic model for the IL-7 and IL-7 receptor system. 



where k+ and k- are, respectively, the binding and unbinding 
rates of the IL-7:IL-7R receptor-ligand system. We assume the 
timescales for this reaction are faster than the timescales for 
changes in receptor numbers, such that we can consider these reac- 
tions to be in equilibrium. This ODE has a unique stable steady 
state: 



R T I 



:fiR T 



The Appendix for Ref. (41) provides estimates for k+ and /c_, 
from which we set /c+ = 1 nMmin -1 and k- = 0.1 min -1 . It is 
reported IL-7 has a molecular mass of around 17 kDa, from which 
we estimate the ratio k-lk+ ~ 1.7ngml _1 . 

A.4. EARLY INTERNALIZATION EVENTS 

In Ref. (40), surface receptor expression was assessed in human 
thymocytes. It is reported, cells in 50ngml _1 IL-7 culture down- 
regulated IL-7R expression. A 20% reduction was observed after 
10 min. Using the above estimate for /c_//c+, we find//l/=5o ~ 0.97. 
Based on this, we can neglect internalization of the unbound recep- 
tor. In the first 10 min, we shall also neglect recycling and inhibition 
of receptor transcription. We assume surface receptor expression 
loss is modeled by the ODE 



drai (t) 
dt 



= 0 - o B mi (t) 



Given initial surface expression levels equal to mi (0), this ODE 
has solution 



mi (f) 



= — + \mi (0) - —1 e ffBt . (A8) 



We assume the previous estimates obtained from the reduced 
model for m\ (0) and 0. That is, we let m\ (0) = 3.6 x 10 4 receptors 
and 0= 1.2 x 10 3 receptors h -1 . Then using the above expression 



for m\{t) equation (A8), we obtain an estimate for cr#. We find 
a B & 1.4 h" 1 . The authors of Ref. (40) estimate the half-life of the 
IL-7 receptor in cells treated with CHX, cultured in 50 ng ml -1 , to 
be approximately 3 h. The expected time to lysosomal degradation 
for a surface receptor is given by 



T2 



_ [cry (1 - ft) + (ju + ju)] gg + gg) + craft (ju + ju) 

CTU (l " ft) ?[/ (§5 + + CTBft (£u + Cu) $B 

^ Uh-^^ + ^ + l^h- 2 

^ 4.6 x 10 3 h" 2 (S B + ( B ) + 1.7 h" 2 ^ ' 



Assuming exponential decay, with a half-life of 3 h, we find 
§£~0.2f£ + 0.3h -1 . Again, without a direct measurement, let us 
set §5=1 h _1 , to obtain an estimate for f b ~ 3.5 h _1 . 



A.5. REMAINING PARAMETERS 

Consider the observation of an approximately 98% reduction 
in surface receptor expression following overnight culture in 
6ngml" 1 IL-7 made in Ref. (39). We let m*\ I=6 = 0.02 m*| /=0 = 



763 receptors. We set the derivatives equal to zero 



0.02 

in the system of ODEs equations (A1-A4), and manipulate the 
resulting system of equations to obtain 



<P<*Bft 



> exp I — 



my 



GU (l -fl)+0 B fl 



<7u(l-fl)%U 0 B fl%B 



(A9) 



exp I —1.9 x 10 receptors 



-) 



0.54, 



3.2 x 10 3 receptors l . 
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Table A1 | Parameter estimates obtained from the mean field model of 
IL-7 receptor dynamics. 



Parameter 


Value 


Units 


0 


1.2 x 10 3 


rech 1 


K 


10 3 


sig 




1 


hn 1 


Hb 


1 


h" 1 




0.14 


h" 1 


OB 


1.4 


h" 1 


kjk+ 


1.7 


ng ml -1 




0.29 


h" 1 


Kb 


3.5 


h" 1 


<P 


0.61 


sig rec -1 hr 1 


X 


0.19 


h" 1 



We let sig be the units of the signaling and rec be the units of surface receptor 
numbers. 



From the steady solutions of ODEs equations (A3) and (A4), 
we find 

(bOnfj 0 W 

ml = raf « 1.9 x 10 z receptors — . 

We use this expression to rearrange equation (A9) in terms of 
m|. This gives 

exp(-^)*9.1^. 

Solving the above expression numerically, we find m|//c ~ 0.1. 
We estimate a value for / based on the observation that following 
culture in 6 ngml -1 IL-7, mRNA levels took approximately 12 h 
to return to 99% of the control levels. The transcription rate is 
given by 0exp(— m±(t)lic). We again assume the IL-7 induced sig- 
nal decays according to the equation m±(t) = m4(0)e~ x \ where 



m4(0) =0.1. Combing these assumptions with the experimental 
observations, we have 

[ m 4 (0)exp(-12 X ) ] 
0 exp = 0.990, 

from which we find x ~ 0.19 h _1 =^ ^ ~ 6.1 x 10 -4 receptors -1 h -1 
k. The parameter k was chosen to be 1000. This choice was made 
from the stochastic model: k = 1000 is the minimum value (to the 
nearest power of 10) such that fluctuations in the signaling quan- 
tity are greater than zero for low ( 10~ 2 ng ml -1 ) concentrations of 
IL-7. Using this value for k, we find 0 ~ 0.61 h _1 . The parameter 
estimates are summarized in Table Al. 

A.5. 1. Changes in the concentration of IL-7 

The functional form 



is used to approximate the numerical solution of m| as a function 
of the concentration of IL-7. Using this function we find a ~ 600 
signaling units and b ~ 0.025 ng ml -1 . In a similar manner we use 
the functional form 

R = c+^— (All) 
e + I 

to approximate the numerical solution of m* as a function of the 
concentration of IL-7. We find c ~ 600 receptors, d ~ 1000 recep- 
tors ngml -1 and e ~ 0.03 ngml -1 . The number of IL-7 molecules 
internalized by each T cell per day is assumed to be the same as the 
number of internalized IL-7:IL-7R complexes per day. We let the 
complex internalization rate, as a function of I, be given by 

24^/jm* (I) « (3 + —i ) 

J 1 2 + 1 V 3 x 10" 2 + l) 

x 10 3 molecules cell -1 day -1 . (A12) 
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